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ABSTRACT 


Passive thermography has been shown to be an effective method for in situ and real time nondestructive evaluation 
(NDE) to measure damage growth in a composite structure during cyclic loading. The heat generation by 
subsurface flaw results in a measurable thermal profile at the surface. This paper models the heat generation 
as a planar subsurface source and calculates the resultant temperature profile at the surface using a three 
dimensional quadrupole. The results of the model are compared to finite difference simulations of the same 
planar sources. 
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1. INTRODUCTION 


Active thermography has been shown to be a very effective method for rapid inspection of large carbon fiber 
reinforced polymer (CFRP) composite structures and is finding increased acceptance. Active thermography uses 
an external heat source to induce a temperature gradient in the material or structure and measures the variation 
in thermal response caused by a flaw. Passive thermography has been shown to be a viable inspection technique 
for hidden defects or damage in structures such as the road or bridge pavement. These thermographic inspections 
rely on thermal gradients in the material or structure, which are due to conditions (such as daytime heating) to 
detect the damage. 


For monitoring the damage growth during carbon fiber reinforced polymer fiber (CFRP) composite fatigue 
testing, passive thermography is almost an ideal technique. Capturing the progression of damage only requires 
an infrared imager monitoring the surface of the composite being tested. The CFRP composite also normally 
has a high enough surface emissivity not to require surface preparation. During the fatigue testing there are at 
least two sources of heating. The first is a result of the dilation deformation of the material and is commonly 
referred to as the thermoelastic effect. For small cyclic strains with sufficiently high frequencies, there is no heat 
transfer to surrounding environment and the mechanism is adiabatic and the time average of the temperature is 
equal to the ambient temperature. It has been shown that it is possible to estimate the strain in the material or 
structure by measuring the cyclic temperature. 


A more interesting heat source during fatigue testing is a flaw in the material. The flaw changes the amplitude 
of the thermoelastic heat generated by changing the load distribution in the vicinity of the flaw. The amplitude 
change is not typically accompanied by a significant phase change. However, if the cyclic loading results two 
surfaces rubbing together, friction generates heat at the flaw. There is phase difference between the cyclic heat 
generation and the cyclic temperature change at the surface, which is a function of the frequency and flaw depth. 
This paper focuses on simulating the heat diffusion from the internal source to the surface. 


While it is desirable to simulate the passive thermography technique, there are few analytical or series solutions 
for the diffusion of heat in a solid and they are only for simple geometries.'_ Numerical techniques such as 
finite difference and finite element methods are well suited for such problems,”*? but can be computationally 
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intensive. Both commercia and noncommercia simulation packages have been successfully applied to 
simulating flash thermography nondestructive evaluation techniques. However, there have been limited attempts 
to simulate passive thermography in composites, in part since passive thermography is primarily a laboratory 
technique. 


A common method for solving the heat equation is to solve the Laplace transform of the heat equation,! then 
invert the Laplace transform to produce a time domain response. There are a limited number of cases where it 
is possible to analytically invert the Laplace transform, however, it is a very powerful technique for one dimen- 
sional multilayered materials when an analytic solution exists for the Laplace transform that can be numerically 
inverted. This is often much faster and more accurate than finite difference or finite element techniques for the 
same configuration, since it is possible to solve for only the times of interest. 


Often this is referred to as the thermal quadrupole method, and is applicable for modeling responses in materials 
and structures in general!® as well as flash thermography. For flash thermography, it has been used extensively for 
simulating the thermal response of multilayer systems with and without contact resistances at the interfaces.!! 1° 
It has also been shown to be applicable for three dimensional configurations, in particular for delaminations 
in composites!’!° by representing the spatial variation in temperature as a cosine series and solving for the 
coefficients. 


This paper follows a similar approach, except the source is interior to the solid and the surface thermal response 
is calculated. The method is illustrated by solving for a temperature response at the surface for a given internal 
source and comparing that to a finite difference solution for the same configuration. While comparisons to 
measurements is preferred, the exact nature of internal sources is difficult to characterize. 


2. QUADRUPOLE METHOD OF SIMULATING THERMAL RESPONSE OF 
LAYERS 


A typical application of the quadrupole method is solving the one dimensional heat equation of multilayer 
systems. For one dimensional problems, a matrix is used to represent the relationship between the temperature 
and flux of one surface of a layer to the temperature and flux at another surface (see Fig. 1). If two of the 
quantities are known (typically the fluxes at the surfaces), then it is possible to solve for the other two. For 
simulating passive thermography, the heat source is not external to the material, but is within the material. The 
next subsection illustrates the quadrupole method in one dimension assuming a flux source at an inner surface. 
This facilitates understanding of the extension of the technique to higher dimensions. 


2.1 Quadrature Method in One Dimension 
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Figure 1. One dimensional configuration for two layer system with a flux at the surface between the two layers. 


A one dimensional solution for the Laplace transform of temperature in a homogeneous material is! 


v(z, 8) = v(0, 8) cosh (<2) F(0, ee (1) 


where K and « are the thermal conductivity and diffusivity, respectively, v(0,s) and f(0,s) are the Laplace 
transform of temperature and flux at z = 0, s is the complex frequency variable and z is the vertical location. A 
similar expression for the Laplace transform of the flux as a function of position is 


f(z,8) = f(0,s) cosh («/*) — v(0,s) sinh («/2) K fe. (2) 


A simple way to express both of the equations is with the matrix formula 
cosh (zq) aed) | v(0, s) | = | v(z, 8) | 3) 
—Kqsinh(zq) cosh (zq) f(0, s) f(z,8) J}? 


where g = Jz . From this matrix equation, the Laplace transform of the temperature at the two surfaces can 
be solved for in terms of the flux, and is found to be 


f (0, s) coth(dq) — f(d, s)csch(dq) 
Kkq 


v(0, 8) = (4) 


and 
csch(dq) f(0, 8) — f(d, s) cosh(dq)) 
- 7 , (5) 
qd 
where d is the thickness of the layer. For flash heating at the front surface (z = 0) and an insulated back surface 


(z = d), the thermal response of the layer can be found by setting f(0,s) = fo, which represents impulse heating 
and f(d,s) = 0, which also can be analytically inverted to give a well known series solution.* 


v(d, s) 


The advantage of this representation is obvious when solving for the thermal response of a multiple layer system 
where each layer is represented by a matrix similar to that given in Eq.3. The matrix is then a transfer matrix, 
and the combined response of multiple layers is obtained by matrix multiplication. For the system shown in Fig. 
1, the two layers are separated into two sets of equations given by 


| cosh (qidi)  =S8bladd | v(0, 8) _ | v(di, 8) (6) 


—Kqsinh (qd1) cosh (qd) f(0,s) f* (d,s) 
and 
cosh (qd2) Se | v(dy, s) | = | u(d, + do, s) (7) 
—Kqsinh(qd2) cosh (qdz) f-(&,8) f(di + da, s) }’ 


where d; and dz are the thicknesses for the first and second layers, f*(di,s) and f~(d1,s) are the fluxes from 
the interface into the upper and lower layers respectively and g = \/s/K. The boundary conditions between 
the two layers represented by this matrix equation are the flux discontinuities due to the heat generated at the 
interface and the temperature is continuous. For no heat generated at the interface, ft(di,s) and f~(dj,s) are 
equal (flux is continuous). 


Considering a special case similar to the passive thermography configuration, with the flux at the front and back 
surfaces equal to the convection losses from the surfaces with an ambient temperature of zero and the ft (dj, s) 
and f~(di,s) are set to fi — fs(s) and fi + f(s), respectively, where f; is an unknown corresponding to the 
heat flux across the surface and f,(s) is the Laplace transform of F;(t), the heat source at the interface. Eqs. 6 
and 7, then become 


cosh (qd) sat) v(0, s) _ | v(di,s) 
—Kqsinh (qd1) oan ais | —h v(0, s) | 7 | fi — fs(s) | (8) 


and 


—Kaqsinh (qd) cosh (qdz) fi + fa(s) h v(di + dz,s) |’ 


where h is the heat transfer coefficient. If f,(s) is given, it is possible to solve this set of equations for temperature 
at either the front or back surface. 


cosh (qd) = Sin(ade) || v(dy, 8) |-| v(dy + do, 8) (9) 


If the thermal properties of the two layers are the same, with a thermal conductivity of K and a thermal 
diffusivity of «, then Laplace transform for vu(0, s) is 


ful) X (VEoosh (day/Z) + fesinh (doy/%) a 
DF os 7 
KO 2h, /= cosh ((di + dz) \/Z) + (8 + (4) x) sinh ((di + d2),/2) 
which can be inverted numerically to calculate the temperature as a function of time. Since there are instabilities 
in the numerical inverse Laplace transform of a sinusoidal function, the numerical inversion is performed for an 


impulse function. The time response for this solution is convolved with the time dependence of the source. The 
time response for F(t) = 1 — cos(2nt) for t > 0 is shown in Fig. 2. The two layers are both 0.1 cm thick, 
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Figure 2. The thermal response calculated from the quadarupole and finite difference methods for a F(t) = 1— cos(2zt) 
excitation 0.1 cm below the surface with and without convection. 


with a thermal diffusivity of 0.005 cm?/sec. Included in the figure are responses for no convection loss (h = 0), 
convection loss of h/K= 1.0/cm and the finite difference solution for both cases. There is reasonable agreement 
between the finite difference solution and the quadrupole solutions. 


For no convection and an impulse heat flux, it is possible to analytically invert the v(0,s) (Eq.10) to give a 
series solution for the thermal response of a one dimension layer. After convolving that solution with cos(wt), 
the surface temperature for a flux source with a time dependence of cos(wt) at d; is found to be 


PRK 


Ucos (0, t) = Kl; +d) (si(t) + cos(wt)se + sin(wt)s3) (11) 


where p is the power amplitude of the source and s1, sz and s3 are given by the equations 
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Figure 3. The thermal response calculated from the quadarupole and series solution for a cosine excitation 0.1 cm below 
the surface. 


A comparison of the solutions is shown in Fig. 3. Shown are the results for a two layer system, where the thermal 
properties of the two layers are the same, the thermal diffusivity is 0.005 cm?/sec, p/K is 1 and dj=0.1cm and 
dg=0.1cm. The solutions are shown for both a frequency of 1 and 2 hertz. The series solution is represented as a 
solid line and the quadrupole solution is represented as a dash line. As can be seen from the figure, the solutions 
are nearly identical for both frequencies. 


For t >> (d, + dz)?/(m?K), the s,(t) is approximately zero, s2 and s3 give the amplitude and phase of the long 
term sinusoidal solution. The amplitude for different depths is shown in Fig. 4. The depth steps are each 0.01 cm 
which is approximately a single ply thickness in a composite. As would be expected, the amplitude significantly 
drops as the frequency increases. As can be seen in the figure, there is a dramatic decrease in the amplitude as 
the depth of the source increases. 


The phase change is shown in Fig. 5. The calculated phase goes from —7 to 7, however, to show the continuous 
variation in the phase, when the phase decreases by more that 7, 27 is added. As can be seen for the figure, 
there is almost a linear relationship between the phase and the depth of the source. As is expected from a simple 
model, the rate of change increases as the frequency changes. 
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Figure 4. The amplitude of thermal response at different depths below the surface. 
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Figure 5. The phase of thermal response at different depths below the surface. 


2.2 Quadrupole Method in Two Dimensions 


Assuming a homogeneous layer with finite lateral dimensions, as is shown in Fig. 6, with no heat flow across 
the vertical edges at x = 0 and x = L, a solution to the Laplace transform of the heat equation in the layer can 
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Figure 6. Two dimensional configuration for two layer system represented with delamination that total blocks flux from 
one layer to the other layer. 


be represented by the expression 


BR 


at) (15) 


N 
v(a, 2,8) = AmUm(Z, §) Cos (= 
m=0 
where L is the lateral width of the layer and a,, = 1 ifm = 0 or m = N —1, otherwise a,, = 2 and the cosine 
series coefficient is t(z,s) is given by 


N-1 rmn 
im(z,8) = aco @nd(tn, 2,8) 608 (Ep™) — Leno AnUn(2, 8) cos (524) (16) 
a al 2N —2 = aN —2 


where the temperature is defined at a set of N evenly spaced locations given by 7, = nL/(N — 1). The flux in 
the layer is given by a similar expression 


N-1 


F(2,2,8) = > amfm(2,8) cos (=), (17) 


m=0 


where panes s) is also given by Eq. 16, if f is substituted for v. Eq. 15 is a solution to the Laplace transform of 
the heat equation if the cosine series coefficients for the temperature have a z dependence and is given by 


ss inh (zm 
Lait Oeneal7.03 (18) 
Kedm 
where gm is given by 
2 
am=\— + (=). (19) 
Ke Ke \ L 


This is similar to the one dimensional equation for the z dependence of the Laplace transform (Eq. 1) and the 
z dependence of flux similar to Eq. 2. The matrix equations for the upper and lower layers are similar to Eq. 8 
and Eq. 9. The two matrix equations for each spatial frequency term are 


a (dma) =sinh(dmdi) | Sm(0, 5) _ _Bn(di, 8) (20) 
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cosh (qmdz2) sn de) | _Um(d1, 8) |-| Om (di + de, s) | (21) 
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where fms(s) is the Laplace transform of the cosine series coefficient the source flux F,(x). There is an assumption 
in this equation that there is no discontinuity in the temperature between the upper and lower layers or that there 
is no contact resistance. The surfaces of the delamination are considered to be in contact enough to generate 
friction and contact resistance of air gaps less than 10 microns is small, therefore is a reasonable approximation. 
The Laplace transform of the cosine coefficient is determined from the equations to be 


jin am ( [Ge cosh (doy) + te sinh (ao /2)) 
FS ap ae own (+ aa,/) ++ Ce) ama (+a) 


By substituting Eq. 22 into Eq. 15, it is possible to find v(2,0,s). This can be numerically inverted to give 
the time response for the surface temperature for a spatially varying subsurface impulse source. This result is 
convolved with cos(wt), to find the time response for a sinusoidal source. An example of the amplitude and phase 
for times significantly greater than (d, + dz)?/(a?«) are of shown in Fig. 7. Phase fluctuations of more than 7 
are removed by either adding or subtracting 27 to maintain a continuous phase profile. The simulation was for 
a source 0.1 cm below the surface in a 0.2 cm thick block, with in-plane and surface normal thermal diffusivities 
both set to 0.005 cm?/sec. A 1 hz source extends from 0.5 cm to 1.0 cm. The convection coefficient is set to 
zero. Included in the figure are the results of a finite difference simulation for the same material properties and 
a convection coefficient of zero. There is reasonably good agreement between the finite difference solution and 
the quadrupole solution. The poorest agreement is in the phase for areas greater than 0.25 cm from the end of 
the source, where the signal amplitude is close to zero. The phase is relatively constant from 0.5 cm to 1.0 cm, 
however, the amplitude is only relatively constant from 0.6 cm to 1.0 cm. This indicates the phase is a better 
parameter for characterizing the size of the source. 
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Figure 7. The amplitude and phase of thermal response for a source that extends from 0.5 cm to 1.0 cm, 0.1 cm below 
the surface. Quadrupole solutions are labeled as Quad and finite difference solutions are labeled FD. 


Since the model assumes no contact resistance at the interface, it is simple to include multiple sources at different 
depths. Results are shown in Fig. 8, where the phase and amplitude are plotted for three 2 hz sources, the first 
0.03 cm below the surface from 0 cm to 1 cm, the second from 1 cm to 2 cm, 0.02 cm below the surface and 
the third from 2 cm to 3 cm, 0.01 cm below the surface. The block was still 0.2 cm thick, however, the thermal 


diffusivity in the in-plane direction is 0.05 cm?/sec and the surface normal thermal diffusivity is 0.005 cm?/sec. 
At the horizontal center of the sources, the phase and amplitude are the same as seen in Fig. 4 and Fig. 5 for 
the appropriate depths of the source. The amplitude profile does not as clearly define the edge of the source as 
the phases as was the case illustrated in Fig. 7. It is also notable that the deeper the flaw, the more ” blurring” 
of the edge. 
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Figure 8. The amplitude and phase of thermal response for three different sources at three different locations from 0.01 
cm to 0.03 cm below the surface. 


3. CONCLUSION 


A method is presented for calculating the time dependence of temperature at the surface from a time depen- 
dent subsurface spatially varying flux. The model is compared to a finite difference measurement of the same 
configuration. The agreement between the two techniques is reasonable. 


A series solution for the one dimensional configuration with a subsurface cosine source is also given. The series 
solution is in excellent agreement with the quadrupole method. The series solution provides an approximate 
time for when the solution becomes totally sinusoidal. 


Future effort will focus on the straight forward process of extending the simulation to three dimensions. It is also 
possible to extend the formulation to multiple layers, which have different thermal properties. The simulation 
outputs will also be compared to experimental results. This is not straight forward, since characteristics of 
internal heat sources are difficult to determine independently. 
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